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Abstract 

The manner in which signals propagate through dense granular systems in both space and time 
is not well understood. In order to learn more about this process, we carry out discrete element 
simulations of the system response to excitations where we control the driving frequency and 
wavelength independently Fourier analysis shows that properties of the signal depend strongly 
on the spatial and temporal scales introduced by the perturbation. The features of the response 
provide a test-bed for any continuum theory attempting to predict signal properties. We illustrate 
this connection between micro-scale physics and macro-scale behavior by comparing the system 
response to a simple elastic model with damping. 



PACS numbers: 45.70.-n,46.40.Cd,43.40.+s, 



The issue of stress and energy transport through dense granular matter (DGM) is of 
significant importance to a number of applications ranging from detection of land mines to 
oil exploration. Due to the importance of this problem, a significant amount of research 
has been carried out in order to understand the basic physical mechanisms involved. The 
major focus of recent work on this issue has been on how the application of static forces 
changes the force structure within a material. One approach involves considering typically 
small size granular samples exposed to time-independent, often point-like perturbations. A 
substantial range of models has been proposed, including diffusive pQ, wave-like [2j [3j H] or 
elastic response [5]. More traditional continuum models [6j [7J) commonly assume an elastic 
or elasto-plastic response. These various models are fundamentally at odds with each other, 
since the basic (possibly continuum limit) equations for these different descriptions are of 
parabolic, hyperbolic, or elliptic nature with respect to their spatial variables. Some progress 
in connecting discrete and continuum descriptions has been reached by realizing that the 
system's response may change depending on the scale or the state of the system: one can 
see wave-like response on short (meso) scales, but elastic response on longer ones [5]. For 
dense systems, theory and experiment [HI E] suggest that an elastic description is best, but 
near jamming threshold, a hyperbolic description may apply [3]. 

Understanding signal (such as a compression wave) propagation through a granular sys- 
tem builds on the approaches used for studies of the static force response but adds the 
additional feature of time-dependence. Within continuum theory, the issue of 'sound' prop- 
agation is often considered via an effective medium approach [7J |9]. A typical model in this 
class explores the response to a spatially independent perturbation which is large compared 
to a particle size, such as the response to a moving piston. Other work has used exper- 
iments and simulations to explore the pressure dependence of the sounds speed, and the 
role of microsctructure including force chains on signal propagation [TT1 [T2l [T3l [T3] • Very 
different issues arise for ID particle chains, which are characterized by nonlinear high-order 
wave-like continuum models [15]. An extension of these results to higher dimensions and 
realistic granular media remains to be carried out. 

This letter concentrates on the response of DGM to space-time dependent perturbations, 
with the goal of gaining further insight into the mechanisms involved in stress and energy 
transport. To the best of our knowledge, this approach has not been considered so far. 
We concentrate on the regime where the imposed frequencies are low, and the wavelengths 



are large compared to the particle size. We then compare the information extracted from 
discrete element (DEM) simulations to expectations based on a simple continuum picture, 
that includes elastic and diffusive behavior. 

We choose a relatively simple granular geometry in two spatial dimensions with the 
granular particles constrained between two rough walls (up-down) with periodic boundary 
conditions (left-right). The walls' position prescribes the volume fraction, v = 0.9. The 
particles are polydisperse disks, with the radii varying randomly in the range of ±5% about 
the mean. For simplicity, we put gravity to zero. The particle-particle and particle-wall 
interactions are modeled using a soft-sphere model that includes damping, dynamic friction 
and rotational degrees of freedom, as explained elsewhere (e.g. [16]). As appropriate for 
2D disks considered here, we use linear springs [T7j. The walls in the simulations are made 
of particles that are rigidly attached, thus creating an impenetrable boundary. The wall 
particles are strongly inelastic and frictional. More detailed explorations of the variation of 
the DEM model parameters, or of the force model itself (e.g., presence of static friction), are 
considered elsewhere [18] . We note that our preliminary results show only weak dependence 
of the response on the details of the force model, or on the parameters. 

The simulations are prepared by very slow compression until the required v is reached. 
After this initial stage, the system is relaxed, the upper boundary is fixed and the lower 
boundary is perturbed as z(x) = z + Asm(ut)sin(kx), where A, u — 2nf, and k = 2n/X 
are the amplitude, angular frequency, and the wavenumber of the imposed perturbation, 
respectively. 

The simulation parameters are as follows. Particle properties: the force constant 
k n = 4000 mg/d, where g is the acceleration of gravity, and m,d are the average mass 
and diameter of a particle; the system particles have the coefficient of restitution e n = 0.9 
and the coefficient of friction fi s = 0.1; the values for the (monodisperse) wall particles are 
e n = 0.1 and fi s = 0.9. System properties: 40,000 particles, with the x dimension being 
L = 250 d. One reason for using such a large system is to reduce the boundary and finite 
size effects. In addition, as we will see below, with the parameters used, some important 
features of propagation are only visible in such a large system. Perturbation properties: 
Amplitude A = 0.6c?; A>d and / <C c*/d (where c* is the speed of sound in the solid) are 
varied. 

Figure [T] shows a snapshot of the simulation domain just after the perturbation of the 
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FIG. 1: Snapshot of the simulation domain soon after perturbation of the lower boundary has been 
activated. The figure shows the forces (normalized by the mean) that the particle experience at a 
given time. The upper wall is static, and the lower one performs standing-wave type of motion. 
40, 000 particles. 

lower boundary has been activated. The color scheme shows the forces on particles, with 
blue (dark) corresponding to low, and green-yellow/light to large forces. We note the force 
chains in the interior of the domain, as typically observed for DGM. 

While the evolution, dynamics, and distribution of force chains is of significant interest 
in understanding signal propagation in DGM [131 EH]; m what follows we consider system 
quantities which can be averaged over a volume which is small compared to the system 
size, but which still includes a relatively large number of particles. Furthermore, we also 
carry out time averaging, choosing an averaging period that is long compared to the particle 
collision time, but short compared to 1//. This averaging procedure is used for calculating 
the quantities which are later compared to the results of a simple continuum model. The 
presented results do not depend on the details of the averaging procedure, as long as the 
general conditions specified above are followed. 

In this work, we concentrate on the elastic energy, which is much larger than the kinetic 
one, to illustrate properties of the propagating signal. Figure [2] shows the elastic energy of 
the granular particles at several different phases during one period of the boundary oscilla- 
tion [19]. We note clear and well-defined wave forms propagating from the lower towards 
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FIG. 2: Energy fluctuations for a DEM simulation where a spatially and temporally harmonic 
perturbation (A = 250<i and / = 30 Hz) is applied at the lower boundary (yellow/light corresponds 
to high energy and blue/dark to low energy). The data in this figure are obtained by subdividing 
the domain in 64x64 cells and than calculating the cell-averages. 

the upper boundary. Analysis of the signal properties such as shown in this figure is one of 
the main points of this work. 

Figure [3] shows the Fourier transforms (FT's), E(z\ f, A), of the elastic energy, E, and of 
the elastic temperature denned as T =< E > 2 — < E 2 > associated with the signal [20]. We 
choose here to present the dominant Fourier mode; the spectral power results are similar. 
These FT's are carried out by expanding the energies accumulated during a given time inter- 
val in the x direction, and then carrying out a time average over a large number (hundreds 
or thousands) of periods of the boundary motion. Time averaging serves to increase the 
signal to noise ratio, and also to ensure that we are not seeing transient effects. The initial 
results, obtained for the first few periods of the oscillations, are discarded. We have further 
verified the steady nature of the results by carrying out selected simulations for much longer 
times. 

Figure [3] shows a well defined signal propagating in the ^-direction. Clearly, E and T, 
which measures local deviations of E from the mean, follow each other closely. Additional 
simulations show that the system-selected wavelengths in the ^-direction are independent 
of the system size [18J. Next we vary / and A of the perturbation, and discuss how the 
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FIG. 3: Dominant Fourier mode (the one imposed by the boundary motion) of the elastic energy 
and temperature as a function of z, the distance from the oscillating wall. 

spectrum of the resulting signal changes. We expect that results such as these will be very 
useful for comparison with any desired continuum model. 

Figure [5] shows E(z; f, A) as A and / are varied. In Figs. [4^,-b we see that a decrease of 
A leads to a significant dispersion, that is, the wave-like property of the signal disappears, 
and the information about the length-scale introduced by the perturbation tends to be lost. 
Note that although A is decreased, it is still large compared to the particle size. Therefore, 
we are not in the regime where finite particle size should be important. Figures [4]>d show 
that an increase of / also leads to the loss of the wave-like signal properties. For a smaller 
/, we see a signal which is weaker and also characterized by a larger z-direction wavelength 
compared to Fig. [3] Therefore, there is just a narrow regime of /'s and A's where a well 
defined signal propagates. For even smaller /'s, the energy input is not sufficiently large to 
be traced in an accurate manner. We note in passing that simulations with an effectively 
infinite A were also carried out, with the results regarding e.g., the speed of propagation, 
consistent with the ones in literature [H] . 

We next compare these results to a model which includes elastic behavior and diffusive 
damping. We chose such a model since static force transmission well above jamming densities 
is, to date, best described in terms of an elastic picture [HI [8], and because we anticipate 
presence of dissipative processes. Assume that the space-time properties of a considered 
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FIG. 4: Dominant Fourier mode of the elastic energy as the frequency and wavelength of the 
perturbation are modified. 

variable E (but similarly for pressure, temperature, or a component of the stress tensor) can 
be described by 



V 2 E - 
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c 2 dt 2 D dt " ' ^ 
which in the limit D — > oo reduces to a linear wave equation often used to describe wave 

propagation in elastic solids, and in the limit c — » oo to a diffusion equation, as considered, 

e.g. in [12] . although for typically much larger /'s. We emphasize here that this model in 

its present form is not meant as a complete description of the simulation results presented 

so far, but instead as a basis for comparison of signal properties that we understand, and 

those what we do not. 

Eq. Q is given in a nondimensional form obtained by choosing L as a length scale, 
and 1/up as a time scale, where u p is some typical imposed frequency (we use u p = 2irf p , 
f p = 30 Hz). We take the diffusion D, and the speed of propagation, c, as constants. Further 
assuming plane wave solution as a zeroth order approximation to more complex waves that 
can be expected in DGM [TH [H)J [2T] in the form E(x, z,t) = E e lu)t e lkx e igz , one obtains the 
following dispersion relation 
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where X = (uj/c) 2 — k 2 . We first discuss general features of such a mode in light of the 
simulation results, and note that the results outlined below apply for a wide range of (con- 



stant) c's and D's. The following predictions are easily verified: (i) for a fixed frequency 
of perturbation /, an increase of k leads to an increase of the dominant wavelength of the 
propagating signal, in agreement with Figs. |3]and|4^; (ii) still for a fixed /, an increase of k 
also leads to larger Im(q) showing that stronger attenuation is expected for shorter A's, in 
agreement with Figs. [3| and |4|d; (iii) for a fixed k, as / is decreased, one expects longer 
emerging wavelengths, in agreement with Fig. We note here that, while the attenuation 
of a propagating signal as a function of driving frequency has been considered before [12j E] , 
we are not aware of any results discussing the dependence of attenuation on the spatial scales 
introduced by the perturbation. One feature of the DEM results which is not explained well 
by the model is the fact that Eq. Q predicts essentially constant attenuation for larger /'s, 
while in Fig. |4ji we see stronger attenuation than e.g. in Fig. [3j consistent with the previous 
work [121 E] ■ One explanation for this difference is that the model predicts shorter emerging 
wavelengths for these high frequencies. These shorter wavelengths may become comparable 
to the particle size, where a continuum model is not expected to apply. 

After finding reasonable agreement between the model, Eq. (jl|, and the simulations, we 
next ask whether the values of D and c deduced by comparison to the DEM data are in 
qualitative agreement with the commonly used ones. For this purpose, we extract the value of 
q from Fig. [3j and using the dispersion relation, Eq. ([2]), obtain c ~ 0.025, D rs 0.01. Similar 
values of c and D can be extracted from the results shown in Fig. [4] and from additional 
simulations using other values of / and A (not shown), typically with the spread of obtained 
values of D larger than the one for c. The value of c can now be compared with the speed 
of sound resulting from elasticity theory, c* = \J E y / ((1 — a 2 )p)/(Luj p ), where E y ,a,p are 
the Young modulus, Poisson ratio and density of the material. These material parameters 
can be extracted from the DEM force model [16] (the force constant there corresponds to 
photoelastic disks such as those used e.g. in [8]), giving c/c* ~ 0.1, in general agreement 
with other works |13j . 

An interpretation of D is more complicated. An estimate based on the particle size 
and some typical shear rate, such as average velocity gradient, underestimates significantly 
the predicted value of D, suggesting that a different mechanism is in place. Alternatively, 
we recall the estimate D « v e l/3 where v e is the velocity of energy transport, and I is the 
transport mean free path, measuring the distance traveled before the direction of propagation 
is randomized; for further discussion regarding applicability of this concept to dense granular 



materials, see, e.g., P, [12]. Let us assume that v e ~ c*. The value of D predicted by the 
model then gives I corresponding to 30 — 40d. It will be of interest to analyze whether such 
a long lengthscale may be related to the correlation length introduced by the force-chain 
structure, the issue which has been a subject of considerable discussion [101 [121 DSl Hi- 
lt will be also of interest to explore how / (and therefore D) vary with, e.g., the volume 
fraction, or the imposed frequency. An additional task will be to analyze the influence of 
dimensionality of the system, since it has been suggested that in 3D the correlation length 
of the force network may be much shorter [12]. We note in passing that, in the regime 
considered here, a diffusion model that was successfully applied to the 3D system driven at 
high frequencies [12], could not explain the main features of the DEM results presented in 
Figs. U and [4j 

While we find consistent results between the simulations and the simple continuum model 
encouraging, we note that much more work is needed. Regarding simulations, we have 
considered propagation for a given volume fraction, therefore not considering explicitly the 
pressure dependence of the speed of propagation. For smaller volume fractions (closer to the 
jamming threshold) the nature of the signal propagation may be modified, and it remains 
to be seen how the continuum model used here would apply in that case. Regarding the 
continuum model itself, although we find that it describes well many features of the DEM 
results, better understanding of the attenuation properties of the signal is needed. The 
degree of agreement with the simulations suggests that for more precise comparison one 
may need to consider frequency-dependent D. The question of coupling of different spacial 
and temporal scales needs to be considered as well. 

We hope that probing DGM with both space and time dependent perturbations, as done 
here, will be utilized in future theoretical and particularly experimental efforts to build a 
more complete picture regarding stress and energy transport in dense granular matter. 
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